Wall influence on dynamics of a microbubble 
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Abstract 

The nonlinear dynamic behaviour of microscopic bubbles near a wall is 
investigated. The Keller-Miksis-Parlitz equation is adopted, but modified 
£SJ to account for the presence of the wall. This base model describes the 

time evolution of the bubble surface, which is assumed to remain spheri- 
cal, and accounts for the effect of acoustic radiation losses owing to liquid 
compressibility in the momentum conservation. Two situations are con- 
sidered: the base case of an isolated bubble in an unbounded medium; 
and a bubble near a solid wall. In the latter case, the wall influence is 
modeled by including a symmetrically oscillating image bubble. The bub- 
ble dynamics is traced using a numerical solution of the model equation. 
^5 Subsequently, Floquet theory is used to accurately detect the bifurcation 

O point where bubble oscillations stop following the driving ultrasound fre- 

* ^ quency and undergo period-changing bifurcations. Of particular interest 

is the detection of the subcritical period tripling and quadrupling transi- 
tion. The parametric bifurcation maps are obtained as functions of non- 
Qh dimensional parameters representing the bubble radius, the frequency and 

pressure amplitude of the driving ultrasound field and the distance from 
I the wall. It is shown that the presence of the wall generally stabilises the 

^ bubble dynamics, so that much larger values of the pressure amplitude 

are needed to generate nonlinear responses. 

00 
(N 

in 1 Introduction 

C*^ When driven by the oscillating pressure field of ultrasound, microbubbles strongly 

CN scatter the incident acoustic waves, and can resonate or fragment [5J [2] • Intra- 

veneously injected microbubbles have been used in clinical practice for over two 
decades [H] because their physical response makes the blood that transports 
J^~j them stand out in the scan relative to the surrounding tissue: this is "contrast" 

rS to the clinician. Modern microbubble contrast agents are coated in an elastic 

shell that retards its dissolution [5J OS] and also affects its properties H7j . 
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Targeted ultrasound contrast agents arc microbubbles in which the shells 
are coated with molecules, usually antibodies, that adhere to specific disease 
markers [IT]. So far, they are not in clinical use. Despite significant research 
into the nature of the signal when microbubbles are bound to targets [TS] , 
at present, there is no clinically- accepted way to tell adherent microbubbles 
from free microbubbles in real time [33J H21 H3J . Since the microbubbles have 
a short lifetime, rapid discrimination of those that are attached to target walls 
from those that have not would be an important step towards clinical practice 
[15]. Suggestions have included filtering based on the fact that the speed of 
adherent microbubbles should be zero [551 H3] , although this does not discrim- 
inate between microbubbles in the blood-vessel boundary layer and ones truly 
adhered. Alternatively, the acoustic radiation force could be used to "push" 
microbubbles; only the free microbubbles would move, enabling the difference 
from bound microbubbles to be discerned [34] . 

For real-time detection of adherent agents, a further suggestion is to exploit 
alterations in the linear response frequency [25) that characterise adherent bub- 
bles. If viscous forces acting in the fluid near the wall are neglected, potential- 
flow theory allows the wall to be replaced with an identical "mirror" bubble 
image that is located symmetrically with respect to the wall plane, see Fig. [l] 
and oscillates with the same frequency, amplitude and phase as the original 
bubble 17]. This mirror- pair oscillates in the symmetric coupled-bubble mode 
[301 126| . If the variable measured is the sound pressure in the liquid, the linear 
natural frequency of the bubble is shifted downwards by a factor of v2/3 ps 0.82 
[IT] . In the medical application that motivated the present work, the walls of 
blood vessels are unlikely to be rigid; in fact the mirror-image symmetric mode 
for a rigid-wall case is one extreme. The opposite extreme is when there is a 
free surface when the bubble and its image oscillate in the anti-symmetric mode 

Q7JHE1B5]. 

Recent experiments with high-speed imaging demonstrated a downshift in 
the linear response frequency of actual targeted contrast agents on biological 
target surfaces [50], broadly consistent with suggestions based on theory and 
experiments on microbubbles larger than contrast-agent sized |25j . With cur- 
rent microbubbles, the diagnostic specificity of such techniques would be poor, 
because the microbubbles vary greatly in size and thus resonant freqency; and 
monodispersed microbubbles are still experimental constructs [7j [3J. Thus, a 
suite of further indicators of proximity to target walls would be beneficial. 
Adding further indicators to the techniques already suggested in the literature 
should improve the sensitivity and specificity of microbubble targeting. 

Detailed research into the nonlinear dynamics of a single, isolated microbub- 
ble began with the work of Lauterborn in 1976 [12] and by 1990 detailed bi- 
furcation diagrams were calculated illustrating the complex behaviour possible 
as driving frequency and amplitude were changed [51]. Recently, the effect of 
acoustic coupling with neighbouring bubbles was examined numerically in [4 
using a coupled-oscillator approach. In the present paper, we expand on the 
earlier suggestions for detection of microbubbles on walls based on linear the- 
ory, to a study of the nonlinear dynamics of a microbubble in the vicinity of 
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Figure 1: Schematic view of the problem geometry. 



a solid wall. In particular, we discuss the non-uniqueness of nonlinear bub- 
ble oscillations that occur for the same driving parameters but different initial 
conditions. 

In addition to the neglect of viscosity inherent in the mirror-image theory, 
we also assume that the distance s/2 between the wall and the bubble centre 
remains constant (in reality, Bjerknes forces would cause this distance to vary 
[TH1I32]). The effects of the bubble shell, which differ widely amongst microbub- 
ble types (and differ amongst models of the shell) [U [24] , are also neglected. 
Such assumptions are made to simplify the calculations and may be refined in 
the future. For the present, the immediate aim is simply to determine if poten- 
tially useful differences exist in the nonlinear dynamics between adherent and 
free microubbles. As noted earlier, the ultimate aim is to derive further criteria 
of wall proximity. 



2 Modified Keller-Miksis-Parlitz equation 



The base model employed in the current study is Keller-Miksis-Parlitz (TUJ [H] 
equation modified to account for the presence of an image bubble, see the last 
term in the right-hand side of equation ([I]), 
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The expression Poo(t) — Pq — Pv + asin(wt), where ui = 2wf e , represents the 
pressure in the liquid far from the bubble, and R(t), Ro, /i, p, K, c, a, a and 
f e denote the instantaneous bubble radius, the equilibrium bubble radius, the 
dynamic viscosity of the liquid, the density of liquid, the polytropic exponent 
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of a gas entrapped in the bubble, the speed of sound in the liquid, the sur- 
face tension of a gas/liquid interface, the acoustic pressure amplitude and the 
driving frequency, respectively. The model considered accounts for the decay of 
bubble oscillations due to viscous dissipation and acoustic radiation. Acoustic 
radiation losses are represented by terms involving the (finite) speed of sound 
c in the Keller-Miksis-Parlitz equation. Bubble oscillations can also decay due 
to thermal energy losses, but such damping is neglected in comparison with 
viscous effects [15] . While the speed of sound is finite in the model employed 
here it is assumed to be sufficiently large so that the phase variation in the 
incoming ultrasound field over the distances of the order of the bubble radius 
are neglected. 

In many bubble-acoustic studies, equations are left in dimensional form. 
However, to reduce the total number of the governing parameters, we make 
equation |l| non-dimensional using the equilibrium radius Rq and the inverse 
ultrasound frequency oj" 1 as the length and time scales, respectively, so that the 
equation is rewritten in terms of non-dimensional bubble radius r = R(t) / Rq 
and time r = Lot as 



f'[(l - fir )r + fiR + Sr 2 ] = (fir - 3) 



W + Rf 



2 r 



_ (M + w) ii±(l_^_ 2Sr . 2 (3) 
- (1 + fir) (M + M e sin r) - M e fir cos r , 



where 



fi= R= -%g, W= 

M = , M e = — §W , S = ^ . (4) 

Each of the non-dimensional groups listed above has a straightforward physical 
meaning. Since in practical applications the driving ultrasound frequency is 
usually fixed, parameter fi, which is the ratio of the equilibrium bubble radius 
and the acoustic wavelength, characterises the bubble size. Parameters R and 
W characterise the viscous dissipation and surface tension effects, respectively. 
They can be treated as inverse Reynolds and Weber numbers. Parameter M 
represents elastic properties of the gas and its compressibility, while M e is the 
measure of the external acoustic excitation. Finally, parameter S is effectively 
the inverse distance between the bubble centre and the wall. 



3 Range of parametric values of interest 

Fluids where microbubble acoustics is of primary practical interest are typically 
water-based. Therefore in order to estimate the values of the governing non- 
dimensional parameters we use the following fluid properties corresponding to 
water at 20°C: c = 1484 m/s, fi = 1CT 3 kg m/s, a = 7.25 x 1(T 2 N/m, p = 10 3 
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Figure 2: Variation of non-dimensional parameters W (dash-dotted line), R 
(dashed line) and M (solid line) with the non-dimensional bubble radius SI: 
(a) SI < 0.02, (b) SI > 0.02. 



kg/m 3 , p v — 2330 Pa. We also assume that the gas trapped inside the bubble 
is air at atmospheric pressure Pq = 10 5 Pa and use the value of K = | for the 
polytropic exponent. 

A typical practical range of the microbubble radii is Ro = 0.5 — 20 /im. The 
frequency used in medical ultrasound imaging is around f e — 1 MHz, and the 
driving pressure amplitude is in the range a = 10 5 - 10 6 Pa. In this study 
we consider bubbles that are assumed to preserve their spherical shape. This 
introduces the natural limitation on the value of s > 2Rq. 

For these physical parameters the values of non-dimensional groups M, W, 
R and SI estimated for bubbles of various radii are linked as 

R= ^ W =^ M e =^° (5) 

SI 2 ' SI 3 ' SI 2 ' e SI 2 ' K ' 

where 

R = 4 ^«1.14xl0- 5 , W = 2 ^«2.79xl0- 7 , 

pc z pc A 

M = w 4.43 x 10~ 5 , (6) 

4.5 x 10" 5 < M e0 = < 4.5 x 10~ 4 . 

for 0.002 < SI < 0.085 (see Fig. [2]), and 10 5 Pa < a < 10 6 Pa. For SI > 
0.008 corresponding to R > 2 pm. the magnitude of M exceeds those of W 
and R, see Fig. [5Jb) . This represents the well-known fact that the dynamics 
of larger microbubbles is mostly determined by the gas elasticity, while both 
viscous dissipation and surface tension play secondary roles. However for smaller 
bubbles the value of W increases rapidly (cubically, see ^) and thus bubble 
dynamics is influenced by the surface tension to a greater degree. The role of 
viscosity for microbubbles of all sizes remains relatively small. 

The other two independent non-dimensional parameters are S and M e . The 
value of S varies from (isolated bubble far away from the wall) to 1/2 (bubble 
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near the wall), while typical values of M e range from to around 0.4 for bubbles 
greater than about 5 /*m, but can be of order 10 for bubbles smaller than 2 /im 
in diameter. 



4 Small amplitude oscillations 

When the non-dimensional driving ultrasound pressure amplitude M e is small, 
the bubble oscillates near its equilibrium state so that its instantaneous non- 
dimensional radius is r(t) = 1 + ri(t), where r\(t) <C 1. Linearizing equation 
([3| about r — 1 we obtain 

(1 + OR + S)r'i + (R + OA)ri + Am 
= -M e (ficosr + sinr) = -M e \/l + O 2 sin(r + <p) , (7) 

where A = (3k — 1)W + 3kM and <f> = tan -1 ft. The solution of this linear 
equation is 

ri(i) = e T o (a sin wo t + 6 cos ujqt) + ao sin(r + 4>o) (8) 

if D = (R — £1A) 2 - 4(1 + S)A < (this condition is always satisfied if the 
viscosity of the fluid remains small as is the case in the problem considered). 
It suggests that the characteristic relaxation time To over which the magnitude 
of the transient solution reduces by the factor of e, the natural frequency lu of 
bubble oscillations, the amplitude ao and the phase shift (f>Q — <j> of the forced 
small amplitude bubble oscillations are given by 



UJq 



«0 




(i + s- A)n 2 



where 



Wo 

m n "o = — ^« 0.0047. 13 

3^°, n»n , 3. Mo 

The approximate values above are obtained by retaining only the largest pa- 
rameters (see Fig. [5] and equations ^ and (j6j)) entering the expressions. These 
formulae provide a straightforward means of identifying the dominant physical 
processes taking place in acoustically forced microbubble oscillations. The value 
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Figure 3: Characteristics of small amplitude oscillations as functions of the non- 
dimensional bubble radius fl for a bubble away (S = 0, solid lines) and near the 
wall (S — |, dashed lines). Intersections of the horizontal dashed lines with the 
Wo lines in plot (b) correspond to subharmonic resonance values of listed in 
Table □ 



of f2o determines the bubble size at which the major physical property deter- 
mining the characteristics of bubble oscillations switches from surface tension 
of the gas/liquid interface to elasticity of the entrapped gas. As noted earlier, 
and now quantified by (13), for the considered fluid properties, only very small 
bubbles with radii smaller than about 1 /im would be affected by the surface 
tension. We also conclude that the bubble oscillation energy dissipation rate 
increases and, subsequently, the relaxation time decreases due to the action of 
two physical mechanisms: a viscous dissipation characterised by R and losses 
due to the acoustic radiation of the bubble, which are proportional to /cfiM, 
yet viscous dissipation always dominates the relaxation process. The proximity 
of the wall (e.g. the increasing value of S) leads to an increase in the relaxation 
time i.e. to the preservation of the oscillation energy due to the reflection of 
acoustic waves from the wall back towards the bubble. This decrease in damp- 
ing owing to the presence of the wall under linear theory has been noted before 

IE]. 

As mentioned in Section [3] the natural frequency of small- amplitude bubble 
oscillations at £1 > ilo (this condition is satisfied for almost a complete range of 
bubble sizes used in practice) is predominantly determined by the gas elasticity 
with the surface tension playing a negligible role. The amplitude of oscillations is 
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directly proportional to that of the ultrasound forcing M e and decreases due to 
the restricting influence of the wall as S increases (when the bubble approaches 
the wall). The phase lag 0o — 4> between the forcing and the bubble response is 
mostly due to the effects of fluid viscosity and acoustic radiation. The proximity 
of the wall forces a bubble to follow the external excitation more closely. 

The variations of characteristics of small-amplitude oscillations with the non- 
dimensional bubble radius are presented in Fig. [3] The overall trends are that 
the relaxation time increases and the natural frequency of bubble oscillations 
decreases with the bubble radius. From a computational point of view it is 
important to note that depending on the initial conditions and the bubble equi- 
librium size it could be required to integrate the governing equations over more 
than 300 forcing periods before a statistically steady solution can be established 
for a bubble located away from the wall. For a near-wall bubble the transient 
time can be greater by at least a factor of two. 

One of the most important sources of information regarding the anticipated 
bubble behaviour is the bubble natural frequency diagram shown in Fig. |3][b). 
The value of luq monotonically decreases with the non-dimensional bubble radius 
and becomes unity at 

(see Table [l] for numerical values) that corresponds to the size of the bubble at 
which it resonates strongly with the external forcing. Typical resonance features 
such as a rapid increase of the bubble oscillation amplitude up to the maximum 
value 

f^M e _ 3^M M e 

a 0ma x« Rq - (1 + S)Rq U&J 

and the switch of the phase lag <fio — <fi from f to — \ are observed at this point, 
see Fig. [3^c, d). For O < fij the amplitude of bubble oscillations decreases 
very rapidly (see Fig. [3](c) ) . A bubble effectively stops responding to an incom- 
ing ultrasound wave and becomes "invisible" in acoustic imaging applications. 
Therefore in our further analysis we focus on the nontrivial bubble behaviour 
observed at fi > fix. We also note that the small amplitude oscillations of 
a bubble near the wall are similar to those of a distant bubble, however the 
resonance shifts towards smaller values of (see Table [lj. This frequency 
shift has been also noted in [2|j] and the potential value of a frequency shift in 
identifying targeted microbubbles via a filtering approach has been recognised 
in [35] . Therefore it is expected that the acoustic signature of initially "silent" 
small bubbles will become more pronounced as they approach the wall. This 
feature can potentially be used for estimating the likelihood of bubbles reaching 
the walls of blood vessels in applications such as targeted drug delivery and 
ultrasound imaging. 

The fact that the non-dimensional natural frequency of bubble oscillations 
remains smaller than 1 for a large portion of the practical bubble size range 
suggests that non-linear subharmonic resonances can occur. The intersections 
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Table 1: The values of parameter Q at which subharmonic resonances with 
frequencies listed in the top row are expected to occur for a bubble away (middle 
row) and near the wall (bottom row). 
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of the horizontal dashed lines with the u>o curves in Fig. Mb) define the sub- 
harmonic resonant values of fl which are listed in Table [l| Of course, the 
subharmonic resonances can only exist owing to nonlinear effects manifested 
when the amplitude is no longer small. As will be shown in the subsequent sec- 
tions, the presence of subharmonic resonances defines what type of oscillations 
are observed. Finally, we note that larger distant and near-wall bubbles follow 
the external forcing quite closely so that the phase lag remains close to zero 
away from the main resonance. 



5 Finite amplitude periodic solutions 

5.1 Floquet analysis 

When the forcing amplitude M e increases, the bubble oscillation amplitude 
becomes large and solutions discussed in the previous section are replaced with 
non-linear solutions ro(r). The solutions are not sinusoidal anymore, yet are 
still T-periodic, where T = 2tt is the period of the external forcing. However 
this remains true only up to a certain critical value of the forcing amplitude at 
which the T-periodic solution undergoes a bifurcation to a different state. To 
determine the exact parametric values for the bifurcation point and the nature 
of the bifurcation we look for a solution in the form r(t) = ro(r)+r'(r), where we 
assume that r'(r) <C ro(r) over the period of one forced oscillation. Numerically, 
tq (t) is obtained by solving equation ([3]) with the periodic boundary conditions 
»"o(0) = tq(T) and ro(0) = f a {T) (see Appendix [8] for details of numerical 
implementation) . To investigate the stability of such a solution we consider the 
linearization of equation ^ about ro(r) 



[(1 - Or )r + OR + Srg]i 



"R, , / , 3. 2 „ .. (3«-l)(M + W) „ . x 
3r + — + 4Sr r + Q M - -rg - R r Q + ± % '- + M e smi 
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(16) 
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and apply Floquet analysis [T!5]. In brief, we solve equation (16) with known 
periodic coefficients depending on ro(r) and ro(r) over the interval r € [0, T] 
subject to two sets of linearly independent initial conditions [r^(0), ^(0)] = [1,0] 
and [r' 2 (0), r 2 (0)] = [0, 1]. The values of the obtained solutions at r = T form a 
monodromy matrix Y 



Y = 



r[(T) 
f[(T) 



r' 2 (T) 
r' 2 (T) 



whose (generally complex) eigenvalues a 12 = <J^ 2 + 2°" 



'1,2 



(17) 



\a x 2 \e i61 ' 2 are Flo- 



quet multipliers. According to the Floquet theorem, the solution of equation 



( 16 ) satisfies the following relationship 



[r'(nT),f'(nT)\ = a n [r'{0), r'(0)] = \a\ n e m9 [r'(0), f'(0)] 



(18) 



Therefore the periodic solution ro(r) is unstable if the magnitude of at least one 
of the Floquet multipliers, maxlcri^l, exceeds unity. At the bifurcation point 
we must have maxlui^l = 1- This condition determines neutral stability of a 
periodic solution with respect to infinitesimal disturbances and the bifurcation 
type is determined by the complex value of the Floquet multiplier with the unit 
magnitude. For example, if a = e 2l7r / n so that a n = a 2lTr = 1 then the solution 
of (|16| will be repeated for the first time after n forcing periods T: 



r'(nT) = r'(0) , f'(nT)=r'(0), r'(mT) ^ r'(0) , r'(mT) ^ r'(0) , 

where m < n, meaning that T —¥ nT bifurcation has occurred. In particular, in 
the situation when one of the Floquet multipliers becomes equal to e l7r = — 1 , 
i.e. when n = 2, a period-doubling bifurcation is observed. 



5.2 Medium amplitude oscillations 

The typical behaviour of Floquet multipliers as the forcing amplitude M e in- 
creases is illustrated in Fig. [4] The magnitude of one of the multipliers increases 
above the value of 1 while the second multiplier remains smaller than unity. The 
phase of the larger of the two Floquet multipliers becomes 9 = tt meaning that 
n = = 2 and thus the period-doubling bifurcation occurs at M e = M ecr . 
Note that a bifurcation of Floquet multipliers seen in Fig. [4] is not related to 
the period-doubling bifurcation of the solution of pi); the latter is determined 
only by the fact that one of the multipliers becomes equal to —1 regardless of 
whether a itself undergoes any bifurcations. Thus the parametric location of a 
bifurcation point is accurately determined by Floquet multipliers that change 
continuously with M e with one of them becoming equal to —1 at the critical 
value of M e = M e cr . 

The period-doubling bifurcations are detected in small to medium ampli- 
tude oscillations for all considered values of non-dimensional bubble radius f2. 
Numerical experiments confirm that the transition between T- and 2T-periodic 
solutions occurs at the same value of M e = M e cr regardless of whether M e is 
gradually increased or decreased. Since no hysteresis is observed we conclude 
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Figure 4: Floquet multipliers for medium amplitude bubble oscillations away 
from the wall (S = 0) at CI = 0.0423. Circles and squares correspond to two 
distinct multipliers. Circles and squares correspond to two distinct multipliers. 
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Figure 5: Bubble oscillations away from the wall (S = 0) at CI = 0.0423 (a) 
before (M e = 0.17) and (b) after (M e = 0.18) the period-doubling bifurcation. 
The solid and dashed lines represent r(r) and r(r), respectively. 
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Figure 6: Oscillations of a bubble near the wall (S = \) at ft = 0.0423 (a) 
before (M e = 0.32) and (b) after (M e = 0.33) the period-doubling bifurcation. 
The solid and dashed lines represent r(r) and r(r), respectively. 

that microbubbles undergo supercritical period-doubling bifurcation in all con- 
sidered regimes. 

The examples of bubble oscillations for 51 = 0.0423 corresponding to the 
equilibrium bubble radius Rq — 10 /im before and after the period-doubling bi- 
furcations occurring at M e cr rj 0.17532 for a bubble away from the wall and at 
M ecr ps 0.32621 for a bubble near the wall are shown in Figs [5] and [6j While 
the qualitative behaviour of the distant and near-wall bubbles remains the same 
we note that the proximity of the wall causes a significant increase in the ultra- 
sound pressure required to induce a supercritical period-doubling bifurcation. 
Equivalently, this means that the presence of the wall has a strong stabilising 
influence on the oscillations of bubbles with the equilibrium radius Rq > 10 /im. 
The pressure amplitude that would have to be applied to cause period doubling 
is nearly doubled by proximity to the wall. In clinical terminology, this implies 
that an ultrasound scanner would have to be set to produce twice the "mechan- 
ical index" in order to see this non-linear response at the wall, a very significant 
change. However, we will show in the following sections that for smaller bubbles 
this trend is reversed. 

5.3 Large amplitude oscillations 

Away from the parametric location of the period-doubling bifurcation the mag- 
nitudes of Floquet multipliers computed for small amplitude T-periodic solu- 
tions (that are accurately approximated by Q) for a bubble located far away 
from the wall remain smaller than unity, see Fig. [7j This means that such oscil- 
lations are stable and the bubble can follow external ultrasound excitation for 
an indefinitely long time (see Fig. |8](a)). However, the phase 9 of the Floquet 
multipliers provides very remarkable information. It turns out that asymp- 
totes to the values ±4?. Thus n = = ±3. Even though no bifurcation is 

O (71,2 

detected, since |o%2| < T the Floquet multipliers strongly indicate that a 3T- 
periodic solution may exist. However, in contrast to the 2T-periodic solution, it 
is disjoint from the T-periodic solution in that sense that it cannot be obtained 
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Figure 7: Floquet multipliers for stable T-periodic bubble oscillations away 
from the wall (S = 0) at fi = 0.0423 indicating the presence of 3T-periodic 
oscillations. Circles and squares correspond to two distinct multipliers. 




Figure 8: (a) Small and (b) large amplitude oscillations of a bubble away from 
the wall (S = 0) at fi = 0.0423 and M e = 0.035. The solid and dashed lines 
represent r(r) and r(r), respectively. 
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Figure 9: Floquct multipliers for stable T-periodic bubble oscillations of the 
bubble near the wall (S = i) at O = 0.0423 indicating the presence of 4T- 
periodic oscillations. Circles and squares correspond to two distinct multipliers. 



from it via a continuous variation of physical governing parameters such as Q 
or M e alone. 

A careful numerical investigation reveals that the governing equation ( 3 1 
indeed admits large amplitude oscillation solutions that are not described by ( 8 1 
even when the external forcing is weak. Consistent with the predictions based 
on Floquet analysis, such large amplitude oscillations have a larger period of 3T 
as is confirmed by Fig. |8jb) . In particular, for the bubble of equilibrium radius 
i?o = 10 fim (f2 = 0.0423), the minimum value of the forcing amplitude at which 
3T-periodic solution is still detected numerically is found to be M e « 0.0342 
(corresponding to a ss 135 kPa) which agrees well with the value estimated from 
Fig. 10(a) in @|. 

Similarly, Floquet multipliers computed away from a period-doubling bifur- 
cation for a bubble near the wall (see Fig. |9j, confirm that small amplitude 
T-periodic solutions are stable, see Fig. [To"ta). However in contrast to the dis- 
tant bubble case, the phase 9 of Floquet multipliers asymptotes to the value 
of #i.2 = ±f so that n = = 4 indicating the presence of 4T-periodic solu- 
tions disjoint from the T-periodic oscillations. Such solutions were indeed found 
numerically, see Fig. [To]^b) . Thus we conclude that the period of large ampli- 
tude bubble oscillations depends on the proximity of the bubble to the wall and 
increases from 3T far away from the wall to 4T next to it. Floquet stability 
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Figure 10: (a) Small and (b) large amplitude oscillations of a bubble near the 
wall (S = |) at fi = 0.0423 and M e = 0.023. The solid and dashed lines denote 
r(r) and r(r), respectively. 



analysis (not detailed here) was performed for the 3T- and 4T-periodic large 
amplitude oscillations; in addition, equation ([3| was solved numerically over 
several hundred forcing periods. This confirmed that the 3T- and 4T-periodic 
oscillations are stable and thus can co-exist with small amplitude T-periodic 
oscillations. 

Having established that the small and large amplitude oscillations illustrated 
in Figs [8] and 10 cannot be obtained from each other by a parametric continua- 
tion, we ask a natural question: what defines the type of the observed oscillations 
in practice? Numerical experiments show that it depends not on the magnitude 
of the forcing, but rather on the initial conditions: large amplitude oscillations 
are typically established if a sufficiently large value of f(0) is specified. Physi- 
cally, this corresponds to a pressure impulse that causes microbubble to initially 
contract with a large speed and then relax to large amplitude long-period os- 
cillations. This is of particular practical significance since clinical ultrasound 
scanners do not apply continuous forcing, but employ a series of discrete pulses. 

These numerical results emphasise that initial excitation can have a long- 
term (in fact, permanent and determining) effect on the observed bubble oscil- 
lation patterns. Therefore the results of the analysis of the bubble dynamics 
based on numerically obtained bifurcation diagrams (Poincare maps) that have 
been a popular tool of bubble dynamics analysis [3TJ [THl IH e.g.] have to be 
interpreted carefully. In obtaining such diagrams it is sometimes assumed that 
the initial conditions are "fully forgotten" once the transients have decayed and 
a periodic solution has been established. However this might not be the case 
due to the existence of multiple solutions of a highly non-linear system ^ each 
having its own basin of attraction [21] . Thus the initial conditions should be 
added to the set of parameters characterising the solutions along with those 
given by Q. Failure to do so may lead to an ambiguous interpretation of the 
computed bubble dynamics. 
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Figure 11: Oscillation map for a bubble away from the wall (S = 0). The 
parametric boundaries for period-doubling and tripling transitions are shown 
by circles and stars, respectively. Vertical dashed lines show the estimated 
positions of subharmonic resonances listed in the middle row of Table [T] The 
dotted curves show isocontours a =const. with the values of a ranging (from 
left to right) from 100 to 1000 kPa at 100 kPa increments. 

6 Bubble oscillation maps 
6.1 Bubble far from walls 

Various types of oscillations that can be experienced by a microbubble located 
far away from a wall are summarised in Fig. |11| The line marked by circles is 
the parametric locus of a period-doubling bifurcation: T-periodic oscillations 
are stable below this line and stable 2T-periodic solutions replace them above 
it. The values of M e corresponding to the period-doubling bifurcation valid 
at least up to 4 significant digits were obtained iteratively by systematically 
applying the Floquet analysis implemented as described in the Appendix. 

The 2T-periodic solution boundary has the characteristic minimum at (CI, M e ) 
(0.029, 0.015) below which no bubbles can oscillate with the half-frequency of 
the driving ultrasound wave. Note that the above value of virtually coincides 
with that at which the natural frequency ujq of small amplitude oscillations be- 
comes equal to | (see Fig. J3] and Table [lj). Therefore the minimum of the driving 
pressure amplitude resulting in the transition to 2T periodic solution is related 
to the first subharmonic resonance observed in the system, an observation con- 
sistent with conclusions regarding the role of resonances in the bubble dynamics 
discussed in |21j . The existence of this minimum could also be qualitatively re- 
lated to the amount of energy required to initiate and maintain oscillations for 
various values of f2. It is intuitively clear that the energy required to start oscil- 
lations (which is proportional to the square of the oscillation amplitude and thus 
to the square of the characteristic bubble radius) decreases with the bubble size 
and so does M ecr . However, as evidenced by Fig. [2] for small values of fl, the 
value of M, the parameter characterising acoustic losses, increases rapidly, and 
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so does the bubble energy loss rate (see the decrease in the relaxation time tq in 
Fig. |3ja) that is inversely proportional to the energy loss rate). To compensate 
for this enhancement of the oscillation energy loss due to acoustic radiation, the 
influx of the ultrasound energy has to increase in order to maintain oscillations. 
Thus the value of M e cr must increase for small bubbles. The minimum in 
the M e cr curve thus corresponds to the bubble size that insures the minimum 
overall oscillation energy loss due to viscous dissipation and acoustic radiation. 

Note that since fi 2 M e = — - we conclude that no bubble can produce a 

pc/ 

half-frequency response if 

a < a 2T ~ 1-26 x lO^pc 2 27.8 kPa. 

Therefore, setting a « an in an experiment and sweeping through a range 
of frequencies / e , one can detect the value of frequency jn at which the half- 
frequency response is first heard. This would provide a straightforward estima- 
tion of the size of the bubble away from the wall 

R ~ 0.029-^- « ~ , 

where f 2 r is in MHz and R is in /jm. For example, for the considered frequency 
f e = 1MHz bubbles for which the half-frequency response would occur first 
would have the equilibrium radius Rq ~ 6.9 pm. 



The line marked by stars in Fig. 11 represents the minimum non-dimensional 
amplitude of the ultrasound forcing at which stable large amplitude 3T-periodic 
oscillations were still detected numerically. These solutions, while remaining 
stable, cease to exist in a catastrophic way (resembling a fold bifurcation) when 
the value of M e is gradually decreased below the starred line. The smallest 
values of M e for which 3T-periodic solutions were still detected were determined 
up to 3 significant digits using the parametric continuation procedure detailed 
in the Appendix. 

The resulting parametric boundary for the existence of 3T-periodic solutions 
has a minimum at (fi, M e ) s» (0.040,0.020). Again the parametric location of 
this minimum is close to the value of fi where the bubble natural frequency ojq 
becomes equal to | of the driving frequency, see Fig. |3|b), and thus it is closely 
linked to the occurring subharmonic resonance. A slight difference between 
the fi location of the minimum and the resonance value reported in Table [T] 
is attributed to a nonlinear shift of the natural oscillation frequency occurring 
for finite amplitude oscillations. Based on the value of M e at the minimum we 
conclude that 3T-periodic oscillations can only be maintained for 

a > a 3T w 3.2 x 10~ 5 pc 2 » 70.5 kPA 

and the bubble size estimation formula based on the period tripling frequency 
fsT becomes 

c _ 9.5 

2tt/3T /3T 



17 




t/T t/T 



Figure 12: Combined bubble oscillations away from the wall (S = 0) at CI = 
0.0423 for (a) M e = 0.18 and (b) M e = 0.25. The solid and dashed lines 
represent r(r) and r(r), respectively. 



where f$T is in MHz and i?o is in /xm. However, given that 3T-periodic solutions 
usually require initial pressure pulse to induce, using this correlation in practical 
experiments may be less convenient than that for 2T-periodic oscillations. 

The 2T- and 3T-periodic solution boundaries intersect at (fi, M e ) ~ (0.0313, 0.0626). 
Therefore it is expected that as the driving pressure amplitude increases the os- 
cillations of smaller bubbles to the left of the intersection point will always expe- 
rience period-doubling before 3T-periodic oscillations can be observed. This fact 
was also illustrated by Figs 10(a) and 12(a) in [4]. The situation is much more 
complicated for larger bubbles. For small forcing amplitudes, in the region below 
the line marked by stars in Fig. [TTJ only small amplitude T-periodic oscillations 
can be observed. In the region between the lines marked by stars and circles, 
both stable T and 3T periodic oscillations can exist, and the type of oscillations 
established is determined by the initial conditions. Above the curve marked by 
circles, T-periodic solutions become unstable and cannot be observed. They are 
replaced by stable 2T-periodic oscillations, see Fig. [5] Again, depending on the 
initial conditions the 3T-periodic solutions also can exist in this region as was 
observed in our numerical computations just above the 2T transition boundary. 
However, for larger values of M e , combination solutions involving both 3T and 
2T components are also observed. In some instances such solutions appear as 



5T-periodic oscillations, see Fig. 12 a). Note that the physical parameters for 
which this plot is generated are identical to those chosen for Fig. j5^b) , yet the 
long-term solutions obtained for different initial conditions differ drastically. In 
other cases the strongly non-linear interaction between these large amplitude 
solutions leads to what appears to be aperiodic oscillations, see Fig. [l2|b), fre- 
quently referred to as chaotic behaviour in literature. Thus we can state that 
the parametric region above both curves in Fig. [TT] is where chaotic bubble os- 
cillations can potentially be observed. This region protrudes towards smaller 
forcing amplitudes near the intersection of the two curves. Therefore the curve 
intersection point defines the size of the bubbles that are most likely to enter 
the chaotic regime away from the wall. For f e = 1 MHz it is Rq ~ 7.4 /zm. 
The higher order subharmonic resonances that are expected to exist for a 



18 



0.4 



(a) 

3 










1 


























1- : 5 \ 



0.1 0.2 0.3 0.4 0.5 




0.031 



0.01 0.02 0.03 0.04 0.05 



Figure 13: Oscillation map for a bubble approaching the wall for fl = 0.0423: 
(a) complete diagram, (b) close up for distant bubbles (computational points are 
shown by symbols, the lines represent a spline interpolation of the data) . The 
parametric boundaries for period-doubling, tripling and quadrupling transitions 
are shown by circles, stars and squares, respectively. 



bubble away from the wall near other vertical dashed lines in Fig. 11 did not 
seem to lead to the appearance of nT-periodic (n > 3) oscillations at least for 
the range of parameters in the figure. We chose not to look for such solutions 
at larger values of M e , since they would be far outside the practical range of 
ultrasound amplitudes that corresponds to the region between dotted curves in 

Fig. [nj 

6.2 Bubble close to a wall 

To investigate the effect of the wall proximity on bubble dynamics, we consider 
an oscillation map produced for il = 0.0423 (corresponding to Rq = 10 /im), 
with S varying from (bubble away from the wall) to | (bubble near the wall). 
The results are presented in Fig. [l3j The only effect the wall proximity has on 
the bubble's tendency to undergo a period doubling bifurcation is quantitative: 
an essentially linear increase of the critical value of M e approximately given by 

M e cr = M e cr0 + dS . 

For fi = 0.0423 and other parameters corresponding to this value (R = 6.4 x 



2|) M eec ,. « 0.0175 and 
that the presence of the 



1CT 3 , W = 3.7 x 10" 3 , M = 2.47 x 10~ 2 , see Fig. 
d s» 0.3. It was noted in the comparison of Figs [5] and 6 
wall can delay the transition to a doubly-periodic regime by a factor of almost 
two in forcing amplitude. This is intuitively expected as the wall is modelled 
by introducing an identical bubble image so that when the bubble approaches 
the wall the external ultrasound radiation effectively has to drive two bubbles 
instead of one. 

In contrast to small amplitude solutions, strongly non-linear large amplitude 
oscillations undergo a qualitative change in the vicinity of the wall: the period 
of such oscillations changes from 3T away from the wall to AT when the bubble 
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Figure 14: Oscillation map for a bubble near the wall (S — |). The parametric 
boundaries for transitions to 2T-, 3T-, 4T- and 5T-periodic oscillations are 
shown by circles, stars, squares and diamonds, respectively. Vertical dashed lines 
show the estimated positions of subharmonic resonances listed in the bottom 
row m Table [Q The dotted curves show isocontours a =const. with the values 
of a ranging (from left to right) from 100 to 1000 kPa at 100 kPa increments. 



approaches the wall, see lines marked by stars and squares in Fig. |13| Such a 
switch occurs between S = 0.01 and S = 0.02, see Fig.[l3)[b) i.e. when the bubble 
is between 25 and 50 radial distances away from the wall. Therefore the influence 
of the wall on the large amplitude oscillations is far-reaching. Another important 
observation is that the forcing amplitude M e decreases along the boundary of 
4T-periodic oscillations as the bubble approaches the wall. Therefore it becomes 
easier to induce such oscillations near the wall. Both these facts suggest that the 
appearance of the quarter-frequencies in the bubble acoustic response spectra 
can be used in practice as an indication of the bubble's approach to the wall. 
In Fig. 14 we present the parametric boundaries for transitions to various 



long-period oscillations of a near- wall bubble. The comparison with Fig. 11 
shows that the proximity of the wall leads to significantly more complicated 
bubble dynamics characterised by the appearance of different oscillation modes 
that depend sensitively on the non-dimensional bubble radius ft. The shape 
of the parametric line showing the locations of period-doubling bifurcations re- 
mains similar to that seen in Fig. 11 for a bubble away from the wall, with the 
minimum detected at the first subharmonic resonance (see the bottom row in Ta- 
ble [TJ. However the figure shows that the behaviour of bubbles of different sizes 
as they approach the wall is not the same. Namely, the excitation amplitude 
at which the period doubling is observed for bubbles with a non-dimensional 
radius ft < 0.026 (Rq < 6 fim) is decreasing as they approach the wall, while for 



larger bubbles it increases. A similar conclusion follows from comparing Figs 11 



and [14] regarding the long-period oscillations however with one drastic differ- 
ence. Near-wall bubbles tend to oscillate with various frequencies that depend 
on their size: the larger the size of a bubble, the longer period of oscillations it 
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has. For example, 3T-periodic oscillations are most profound near the second 
subharmonic resonance (stars in Fig. 14). They are replaced with AT- (squares) 
and 5T- (diamonds) periodic oscillations as increases, but such long-period 
oscillations are not expected to be seen for forcing amplitudes commonly used 
in medical practice. As noted earlier, the tendency of a near-wall bubble to un- 
dergo longer period oscillations is apparently linked to the presence of a bubble 
image in the model that effectively increases the bubble inertia; note that the S 
term plays the role of an additional mass in the left-hand side of equation 

It is also instructive to note that irregular oscillations can be observed at the 
values of M e smaller than those shown by squares and diamonds in Fig. |T4| ) for 
some initial conditions. They could be mistaken for "chaotic bubble behaviour" . 
In reality they appear to be just transient solutions that eventually settle to 
simple T-periodic oscillations. However this may take a very long time (up to 
20,000 of forcing periods in our computations). The fact that they eventually 
decay while regular multi-period oscillations persist confirms once again that 
the appearance of the latter is due to subharmonic resonances which facilitate 
absorption of ultrasound energy by the resonating bubbles. 



7 Conclusions 

The behaviour of a microbubble approaching a solid wall and subjected to ul- 
trasound forcing has been investigated using a modified Keller-Miksis-Parlitz 
equation [TUl [H] as a base model. Floquet analysis has been applied to in- 
vestigate the stability of the small amplitude solutions and predict the regions 
of existence of various fully nonlinear large amplitude oscillations. It has been 
shown that microbubble response to acoustic forcing can consist of T-, 2T- and 
3T-periodic oscillations (T is the period of ultrasound forcing) when the bubble 
is away from the wall and T-, 2T-, AT- and 5 T-periodic oscillations when the 
bubble is near it. We also showed that a bubble starts "feeling" the presence of 
the wall and changes its acoustic response when it as far from the wall as 25-50 
radial distances away. It is found that in all cases 2T oscillations appear as a 
result of supercritical period doubling bifurcations of T-periodic solutions, while 
longer-periodic oscillations are completely disjoint from the T- and 2T-periodic 
solutions and usually require a pressure impulse to be initiated. 

Although several physical details were neglected for simplicity, such as the 
presence of a shell encapsulating medical microbubbles, the Bjerknes force that 
could cause a bubble on a wall to deform, and compliance of the wall material, 
the existence of clear differences in the dynamical-systems behaviour owing to 
the presence of a wall suggests that more detailed studies are warranted. The 
results of future investigations could be used to estimate both the size of the 
bubble and its proximity to the wall based on the spectrum of the bubble's 
acoustic signature. 

The demonstrated existence of multiple oscillatory solutions also strongly 
indicates that bifurcation diagrams (Poincare maps) that are frequently used in 
studies of bubble dynamics could be misleading unless initial conditions used to 
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obtain such diagrams are clearly stated and added to the set of the governing 
parameters characterising the problem. 



8 Appendix. Numerical aspects 

The major value of Floquet analysis is in its ability to predict very accurately 
various period-changing bifurcations. However this requires the explicit knowl- 
edge of the periodic solution (V (t), r (r)) stability of which it aims to investi- 
gate. One may try to use a "brutal force" approach by integrating the governing 
equation (|3| starting with some arbitrary initial conditions over a sufficiently 
long time to allow all transients to decay and a periodic (limit cycle) solution 
to establish. However this can only work if this solution is stable. 

Yet to determine the bifurcation point iteratively one needs to compute Flo- 
quet multipliers and thus know the periodic solution (ro(r), ro(r)) in parametric 
regions where this solution becomes unstable and thus cannot be obtained us- 
ing forward time integration. Therefore instead of solving ([3| as an initial value 
problem one needs to view it as a periodic boundary problem with a specified 
period (the forcing period or its integer multiple). Solving such a boundary 
value problem is done iteratively and the convergence of iterations depends sen- 
sitively on the initial guess. For T-periodic oscillations it is readily available 
from a large time limit of ([8| and the iterative determination of one period- 
doubling bifurcation point takes just a few seconds of CPU time. 

However this analytic solution cannot be used as an initial guess for large 
amplitude long period oscillations. In fact, it is not even known whether such 
solutions exist as they are fully disjoint from the T- and 2T-periodic solutions 
naturally linked to Q. Therefore first we used the analysis of the phase 9 of 
Floquet multipliers (see Figs [7] and [9]) of the T-periodic solutions to establish 
parametric ranges where fully nonlinear large amplitude long period solutions 
can possibly exist. Then for a so selected set of governing parameters we tested 
a number of empirically chosen "pressure impulse" initial conditions to obtain 
a stable large amplitude solution using forward time integration. Numerically, 
this was done by specifying a sufficiently large (typically of order 10 _1 or even 
1) negative value of ro. 

Once the first long term nT-periodic (n > 2) solution was obtained in such 
a way, similarly to Q offering a suitable initial guess for solving a periodic 
boundary value problem for small amplitude solutions, it can be used as a 
natural initial guess for a periodic boundary value problem for large amplitude 
oscillations. However we found that solving it for rtT-periodic oscillations still 
encountered significant numerical difficulties due to the cusp-like singularity 



developing near the minima of the bubble contraction curves (see Fig 12). This 
singularity frequently results in the divergence of iterations due to the local 
loss of numerical approximation near the cusp. Since the location of such a 
cusp changes during iterations, adaptively increasing the local grid density near 
it to recover the approximation accuracy becomes a technical challenge. For 
this reason in the current study we chose not to solve a periodic boundary 
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value problem for long period nonlinear oscillations. Instead we used a more 
robust forward time integration e.g. implemented in Matlab's function odel5s 
for numerically stiff problems for all large amplitude oscillations since it has a 
built-in automatic algorithm for reducing a computational step when the time 
derivative becomes large. 

As discussed above, only stable oscillations can be obtained using this ap- 
proach so that an automatic iterative search of the transition point based on 
Floquet multipliers that was successfully implemented for small amplitude oscil- 
lations cannot be guaranteed to work for large amplitude nT-periodic solutions. 
Therefore instead, once the first of such solutions was obtained as discussed 
above, the governing physical parameters were gradually changed and the final 
values of (ro,ro) from a previous forward time integration run were used as 
initial values for a new run in order to trace the variation of large amplitude 
solutions with M e , S and f2. The length of each such run was between 900 and 
4000 forcing periods (between 50 and 300 relaxation time units To) to allow for 
statistically steady oscillations to establish. The longest runs took up to 4 min 
on a standard desktop computer. 

were 



The lines representing nT-periodic solutions in Figs 11 13 and 14 



obtained in this way and thus they correspond to parametric values at which 
the respective long-period numerical solutions cannot be found any more by 
forward time integration using the above parametric continuation procedure 
towards smaller values of M P . 
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